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ABSTRACT 

This paper concerns the interaction between non-axisymmetric inertial waves and 
their corotation resonances in a hydrodynamical disk. Inertial waves are of interest 
because they can localise in resonant cavities circumscribed by Lindblad radii, and as 
a consequence exhibit discrete oscillation frequencies that may be observed. It is often 
hypothesised that these trapped eigenmodes are affiliated with the poorly understood 
QPO phenomenon. We demonstrate that a large class of non-axisymmetric 3D inertial 
waves cannot manifest as trapped normal modes. This class includes any inertial wave 
whose resonant cavity contains a corotation singularity. Instead, these 'singular' modes 
constitute a continuous spectrum and, as an ensemble, are convected with the flow, 
giving rise to shearing waves. Lastly, we present a simple demonstration of how the 
corotation singularity stabilizes three-dimensional perturbations in a slender torus. 
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'. 1 INTRODUCTION 

V.O Inertial waves, also called r-modes (Korycanksy and Pringle 
O ' 1995, Ogilvie 1998) or g-modes (Wagoner 1999, Kato 2001a), 

^ , are one of a rich assortment of oscillations exhibited by 

hydrodynamical models of accretion disks. They may be 
distinguished by a number of interesting properties: first, 
f*"^ . they comprise motions that are effectively incompressible, 
and which strongly couple the vertical and horizontal veloc- 
. ,_| \ ities. Second, their oscillation frequencies take values much 
i less than the local epicyclic frequency, which suggests that 
'_h inertial waves figure prominently in a disk's response to 
^ low frequency forcing (Balbus 2003), such as might issue 
from an embedded protoplanet. Third, these modes are the 
only waves that propagate in the vicinity of their corota- 
tion resonances, i.e. the radii where their pattern speed is 
zero. Finally, in principle they may become trapped between 
their two Lindblad resonances and hence manifest as a dis- 
crete suite of standing modes (see, for example, Okazaki et 
al. 1987, Perez et al. 1997). 

Inertial waves have been invoked to explain various as- 
trophysical phenonomena, but of special interest is their 
potential role in the excitation of the quasi-periodic oscil- 
lations (QPOs) that are observed in the X-ray spectra of 
some black-hole disk systems (e.g., Wagoner 1999; Kato 
2001a, McClintock and Remillard 2003). The idea is that 
since the inertial modes appear as trapped standing waves, 
they (or their resonant interactions) may exhibit oscillation 
frequencies in accord with the observed frequencies. Mean- 
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while, a variety of instabilities have been promoted in order 
to explain the large amplitudes necessary for the detection 
of these frequencies (Abramowicz and Kluzniak 2001, Kato 
2001b, 2003b, Ferreira and Ogilvie 2008, for example). 

Though axisymmetric waves (Okazaki et al. 1987) and 
certain non-axisymmetric waves (Perez et al. 1997) have 
been demonstrated mathematically to form discrete trapped 
standing waves, this behaviour has not been established sat- 
isfactorily for general three-dimensional modes. In this paper 
we show, with both WKBJ techniques and direct methods, 
that in most cases, three-dimensional inertial waves do not 
manifest this way. This is because most non-axisymmetric 
3D modes possess a corotation singularity within their reso- 
nant cavity. While earlier work (Kato 2003a, Li et al. 2003) 
has shown how the corotation singularity damps incident 
wave fluxes, it has not been generally appreciated that the 
corotation singularity also very likely prohibits the forma- 
tion of trapped normal modes themselves. The present paper 
is concerned with investigating the dynamics of this prob- 
lem. In addition, we study the continuous spectrum that 
issues from the corotation singularity, showing its close rela- 
tionship with inertial shearing waves (Johnson and Gammie 
2005, Balbus and Hawley 2006). Lastly, we demonstrate how 
the wave absorption at corotation impedes unstable three- 
dimensional modes in a slender torus. In so doing, we deepen 
our understanding of the complicated and varied linear re- 
sponses available to accretion disks. 

The outline of this paper is as follows. Section 2 presents 
the governing equations of an incompressible disk in a shear- 
ing sheet and summarises the main qualitative features of 
travelling inertial waves near the corotation point. Section 3 
demonstrates the absence of trapped modes around corota- 
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tion by solving the boundary value problem. In the following 
section we provide details of a WKBJ analysis which sup- 
plies the same result in a more realistic semi-global model 
but in the limit of large vertical wavenumber. Section 5 out- 
lines the relationship between shearing waves and the contin- 
uous spectrum induced by the corotation singularity, while 
Section 6 describes instability in a 3D slender torus. Our 
conclusions are drawn in Section 7. 



2 PRELIMINARIES 

Our treatment is restricted to inertial waves whose resonant 
cavities contain the corotation singularity. Such a configura- 
tion is typical for 3D inertial waves in both Newtonian and 
(general) relativistic disks, with only a narrow band of non- 
axisymmetric r-modes in relativistic disks exhibiting a differ- 
ent configuration of resonances (see Perez et al. 1997, Wag- 
oner 1999, Kato 2001a, 2001b). The more common modes 
were first studied in detail by Kato (2002), who suggested 
that the corotation point could give rise to a form of the 
Papaloizou-Pringle instability. Subsequent WKBJ and nu- 
merical analyses, however, demonstrated that an incident in- 
ertial wave would be strongly absorbed at corotation (Kato 
2003a, Li et al. 2003, but see Drury, 1985, for a much earlier 
discussion). These results convinced researchers that such 
trapped 3D modes must decay. As a consequence, the em- 
phasis of much 'diskoseismology' has moved on to other 
mechanisms of instability, typically involving a paramet- 
ric resonance (Abramowicz and Kluzniak 2001, Kato 2003b, 
Ferreira and Ogilvie 2008). However, a central question re- 
mains unaddressed: How is it possible that trapped standing 
waves (growing or decaying) form at all in the vicinity of a 
strong wave absorber, such as the corotation singularity? 

The WKBJ analysis of Kato (2002, 2003a; hereafter 
K02 and K03a), assumes the existence of such trapped 
modes and computes their damping rates. The theory in- 
corporates neither explicit nor implicit dissipation, which al- 
lows a growing mode for each decaying mode present, due to 
time symmetry. But because the problem exhibits a branch 
cut, this symmetry must be broken in any calculation. At 
the outset, one must assume either growth or decay of the 
modes to be studied. Henceforth, care must be taken to as- 
sure self-consistency. Both K02 and K03a assumed that their 
modes grow in their calculation. Their final result, however, 
was that the modes decay. This was treated, not as an in- 
consistency to be rectified, but as indication the modes do, 
in fact, decay. It is the view of the authors that the correct 
interpretation of the inconsistency is that there exists no 
trapped modes at all. The only class of modes that could in 
principle resolve the contradiction are strictly neutral. But 
such modes are singular, possessing undefined radial deriva- 
tives at the critical point. 

In Section 4 we revisit the analysis of K02 and K03a 
and reestablish the above result in some detail. Before do- 
ing so, we demonstrate the main ideas using a simple lo- 
cal model, the shearing sheet, and very straighforward tech- 
niques. In the shearing sheet limit, the trapped wave prob- 
lem can be solved analytically, and hence is more readily 
understood. This limit neglects global phenomena such as 
curvature terms and edge effects, as well as general rela- 



tivity, but it includes all the essential physics: a corotation 
singularity with two Lindblad resonances on either side. 

We begin with the basic equations in the shearing sheet 
limit, and then give a brief summary of travelling inertial 
waves. This will help us interpret the results of the follow- 
ing section in which the absence of normal modes is demon- 
strated. 

2.1 Governing equations 

The shearing sheet (Goldreich and Lynden-Bell 1965) is 
designed to study local disk behavior. We make the ad- 
ditional assumptions that self-gravity is unimportant, and 
that the fluid is incompressible. Also, we work in the 'cylin- 
drical limit': the equilibrium density is unstratified, inde- 
pendent of both cylindrical radius (r) and height (z). It can 
be shown, however, that the vertical stratification plays no 
role in the qualitative behaviour of our problem. For exam- 
ple, an isothermal model yields analogous results, at greater 
mathematical cost (see Appendix A) . 

The governing equations in the shearing sheet read: 

dtu + u- Vu + 2n e* x u= -V$ - -Vp, (1) 

P 

V ■ u = 0, (2) 

where u is velocity, $ is the tidal potential, p is pressure, 
and p is density (assumed a constant). Additionally, fio = 
n(ro) is the local orbital frequency of the disk, where ro is 
the radius at which the shearing sheet is anchored. The x- 
coordinate represents radius (dx = dr) , the y coordinate the 
azimuthal direction (dy = rdff), and the z coordinate the 
vertical direction. The sheet undergoes a background linear 
shear flow 

uo = — <7^o x e y (3) 

where q = —(din Q/d\nr)o. If the disk is Keplerian, q = 3/2. 
Henceforth, the subscript on fi will be suppressed. 

2.2 Linear perturbation equations 

The governing equations admit by construction the homo- 
geneous equilibrium u = uo, p = po, and p = po- Let us 
now introduce a small disturbance so that u = uo + u', and 
p = Po + p' ■ The equations governing the linear evolution of 
such a pertubation are: 

d t u x — qQx d y u x — 2Q,u' y = —d x h', (4) 

dtu' y — qClxdyUy + |(k 2 /Cl)u' x = —d y h', (5) 

d t u' z — qQx d y u' z = —d z h', (6) 

with the incompressibility restriction: 

d x u' x + dyu'y + d z u' z = 0. (7) 

Note that the pertubed enthalpy h' = p' / po takes the place 
of the perturbed pressure p . The epicyclic frequency is de- 
fined through K 2 = 2(2 - q)Q. 2 . 

The perturbations are Fourier decomposed in space so 
that each mode is oc exp(ik y y + ik z z). The adoption of sinu- 
soidal variation in z is permitted by the cylindrical approxi- 
mation. After some manipulation we arrive upon the central 
equation of the linear analysis for a single mode (k y , k z ), 

(dt - iq0.xk y ) 2 (dl - k±)u' x + n 2 k 2 u' x = 0, (8) 
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where k\ = k 2 + k 2 , and u' x now represents a Fourier ampli- 
tude. At this point we perform a temporal Fourier transform 
and assume the equation admits discrete modal solutions, 
i.e. we take u' x oc e -1 "', where uj is a (complex) frequency. 
We set 

uj = uj r + ia, 

for uj r and a real. The sign of a must be specified now so 
as to negotiate the branch point that arises from the invis- 
cid dynamics. Without loss of generality we set a ^ 0, as 
the inviscid problem supports either neutral modes or grow- 
ing/decaying pairs (due to the time symmetry). Equation 
((8| is now 



where 



— 2 o2 / i / 2 1 2 j 2 — 2\ / r\ 
UJ O x U x + (K K z — k±UJ )U X = U, 



uj — uj + qQ x ky . 



(9) 



(10) 



2.3 Wave-shape equation 

In order to achieve a simpler version of our working equation, 
a shifted dimensionless space variable is adopted, 



= k± [ x + 



qQ. ky 



(11) 



Hence we anchor the shearing sheet upon the radius in the 
disk at which the mode is convected by the shear flow. The 
star is dropped hereafter for notational convenience, in ad- 
dition to the subscript x and prime on u' x . The central wave- 
shape equation of the analysis is consequently 



dZu + 



ii 2 



(x- 



0. 



where 



qQ,k y 



(12) 



(13) 



It will be assumed, without loss of generality, that k y > 
and so x c < 0. The dimensionless parameter j3 is defined 
through 



qQ 



(14) 



In fact, p 2 is directly analogous to the Richardson num- 
ber (Booker and Bretherton, 1967, and see Li et al., 2003), 
and measures the ratio of the stabilizing effects of rotation, 
embodied in the term K 2 (k z /k 2 ), against the destabilising 
effects of shear. 

The governing wave-shape equation possesses a singu- 
larity at x — ix c (the corotation point). However, the solu- 
tions are regular for real x provided that x c 7^ 0, meaning 
that the solutions must grow (or decay). Strictly neutral 
modes (for which x c = 0) are undefined because they have 
singular derivatives. With a > 0, and x c < 0, the singularity 
always lies beneath the real axis, which serves as our contour 
of integration. Thus the branch cut extending to the left of 
the singularity always lies beneath the integration contour. 

The singularity corresponds to the corotation reso- 
nance, and arises because the fluid is inviscid; it vanishes 
when dissipation is added. Note also that the singularity 
depends on k z 7^ (from Eq. (|14[) ), and is not present for 



two-dimensional modes. The two Lindblad resonances occur 
at 

x = ±/3 + ix c , 

which are arranged symmetrically around the corotation ra- 
dius. (This symmetry is broken when curvature effects are 
included.) The resonant cavity is defined by \x\ ^ (3. Den- 
sity waves can propagate outside of this region, but because 
the fluid is incompressible, the 'evanescent zones' extend to 
positive and negative infinity. 

Trapped standing waves must decay far from corota- 
tion in the evanescent zones: \u\ — » when \x\ — > 00. Equa- 
tion (I12p is then a one-dimensional eigenvalue problem with 
eigenvalue a. In fact, Eq. (|12|l is a form of the Bessel equa- 
tion and so the general solution can be written in closed 
form immediately 



u = ai \fz 



) + bi \/x — ix c I- V (x — ix c ) (15) 



where ai and 61 are constants, I v (z) is a modified Bessel 
function of first kind, and the order v is defined through 



1 



= -yr^r, 



(16) 



(see, for example, Abramowitz and Stegun, 1972). This is 
particularly advantageous as, unlike other work, we need 
not resort to WKBJ or Frobenius approximations. The for- 
mer corresponds to the limit of large imaginary order v (i.e. 
f3 3> 1) and the latter to small argument x. Two classes of 
solution can be distinguished. When /3 > 1/2 (corresponding 
to modes with greater wavevector pitch k z /k y ) the order of 
the Bessel functions v is imaginary and strongly oscillatory 
solutions are obtained. When /3 ^ 1/2 the profiles resemble 
hyperbolic trigonometric functions. 

Before we attack the boundary value problem we discuss 
free travelling waves within the resonant cavity. This will 
illuminate some of the properties of the standing waves near 
corotation and help us understand in more physical terms 
the results of Section 3. Some of this analysis also appears 
in Booker and Bretherton (1967), Vishniac and Diamond 
(1989), and Li et al. (2003). 



2.4 Free inertial waves near corotation 

We first discuss the properties of wave-like solutions to <(9j 
far from corotation uj = by making the familiar WKBJ 
ansatz, u — R^e'^^^, where ip is a phase function that 
varies rapidly with x in comparison with R. We make the 
identification k x = d x ip and also assume that ky ^ k z ^ k x . 
This ansatz is substituted into © and we collect terms of 
order k x . These assumptions lead to the dispersion relation 
for inertial waves, 



k z 



(17) 



where k 2 = ki 



~\~ ky -\- k. 



(see also Balbus 2003). The group 



velocity of these waves with respect to the Doppler-shifted 
frequency uj is given by 



V = V k S 



UJ 
fc2 



-k y , (k x + ky)/k z 



(18) 



where k = (k x , k y , k z ). Clearly, k ■ V = and the group 
velocity is perpendicular to the phase velocity. Moreover, 
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the sign of the radial component of the group velocity V x 
depends on both k and the sign of uj, and consequently the 
location of the wave relative to the corotation point. Inside 
corotation (uj < 0), the group velocity points in the same 
radial direction as that of the phase velocity. But, outside 
corotation (uj > 0), the group velocity points in the opposite 
radial direction to the phase velocity. 

Next we treat the fate of the inertial waves as they ap- 
proach corotation, assuming that the group velocity calcu- 
lated above does not deviate significantly in this region. Near 
this radius a Frobenius solution may be employed. From 
Eq. Q5) , to leading order the approximation is 

u = ai(x- ix c ) 1/2+l/ +bi(x- ix c ) 1/2 - v , (19) 

where oi, &i and v are the same quantities introduced ear- 
lier. We adopt the WKBJ approximation and thus assume (3 
is large; now v « i/3. Consequently, a local radial wavenum- 
ber can be defined for each of the two solution components 
above. In particular, 

u = \J x — ix c (aie mnlx - x ° l + bie -i0^\*-*c\\ . (20) 

The (x — XcY 13 component possesses k x = fH/x, and con- 
sequently we identify the a\ function as a left-going wave 
(moving from positive to negative x). This is because when 
x > then k x > 0, and so, according to Eq. (|18[) . the ra- 
dial component of the group velocity V x is negative. When 
x < 0, then k x < 0, and V x will also be negative. A sim- 
ilar argument establishes that the b\ function is the wave 
propagating to the right (moving from negative to positive 
x). 

When one of these wave components passes through 
corotation, it will be strongly absorbed, suffering an am- 
plitude damping proportional to e -7 ' 13 . This is clear from 
the analytic continuation of the Frobenius solution Eq. (|19[) 
through the corotation region, noting that the singularity 
x — ix c lies beneath the real line, and thus beneath the 
integration contour. For large /3, i.e. large wavevector pitch 
k z /ky, the incident waves barely penetrate the corotation 
barrier. For general j3, when each wave passes through coro- 
tation from its proper direction, its amplitude will acquire 
a factor 

exp(i in + iirv ), (21) 

from (|19[) . For ^ 1/2, the alteration corresponds to a 
phase shift of -k(v + 1/2). But it nevertheless gives rise to 
a damping of angular momentum and energy flux through 
the corotation point (Li et al. 2003). 

It is this damping that has led researchers to discount 
the possibility of instability in this context: the corotation 
resonance is a powerful barrier against the communication of 
wave motions, plundering energy from potentially growing 
modes. We now show that this barrier not only precludes 
instability, it precludes standing modes entirely. 



3 ABSENCE OF TRAPPED 3D INERTIAL 

WAVES IN THE SHEARING SHEET MODEL 

In this section we take the exact solution (|15p to the gov- 
erning wave-mode equation (|12[) and show that it cannot 
satisfy the two decaying boundary conditions at \x\ — > oo. 



Physically, the reason for this is that the standing mode con- 
sists of a precise balance of two traveling waves (the left and 
right-going WKBJ waves of the previous section) and the 
corotation point will heavily damp one component relative 
to the other. Consequently, the correct balance of compo- 
nents on one side of corotation will never correspond to the 
correct balance on the other side of corotation. 

Mathematically, the result is straightforward. The 
boundary condition at x — > oo forces u to be proportional 
to the Bessel function of second kind K v (x — ix c ) for all x, 
as long as x c 7^ 0. This function, however, diverges when 
x — > —00 and so no solution exists. 

Even though the mathematics is clear, we develop this 
problem in some detail because it sheds light on the physics 
of the mode prohibition. In particular, it shows that the 
problem issues from the fundamental interaction between a 
wave and its corotation singularity, and is thus quite general 
— not an artefact of the shearing sheet model. After all, 
the shearing sheet also prohibits k z = modes (for which 
there is no singularity) because of its simple linear shear (by 
analogy with Case 1960). Moreover, the lengthier treatment 
below helps us deal with the nontrivial x c — case, which 
sets up a treatment of 'singular modes' and the continuous 
spectrum. 

3.1 Boundary conditions 

We apply the boundary condition ii| — > as — > 00 upon 
the general solution to (|12[1 . 

it = ai x — ix c Iv(x — ix c ) + bi yx — ix c I-„(x — ix c ). (22) 

Normally these constraints would lead to two algebraic equa- 
tions for 61/ai and x c . Determining the latter is the 'eigen- 
value problem', or 'quantisation condition'. As we shall see, 
such an equation does not appear in this problem. 

Throughout this section, we take (3 large (WKBJ waves) 
and so make the approximation v « if3. Each solution com- 
ponent above can then be regarded as a travelling wave: 
the Ii@ component corresponds to (x — x c ) 113 and the left- 
travelling wave, while the I-ip component corresponds to 
the right-travelling wave. These two must be balanced ap- 
propriately in order to form a localised standing wave. 
Lastly, recall that k y > and thus x c < 0. The sign of x c is 
important in what follows as it determinines the asymptotic 
forms of the Bessel functions far from corotation. 

Consider first the limit x — > 00. The asymptotic forms 
of the two solution components can be derived using 

where X = x — ix c (Abramowitz and Stegun 1972). To force 
the solution u to decay when x is large and positive, the 
coefficient of exp(x) must be set to zero, and this gives us 
the relation 

bi/01 = -1. (23) 

This combination means that u is proportional to 
K v (x — ix c ), the associated Bessel function of second kind 
(mentioned above). 

Consider next the other side of the critical point in the 
limit x — ► —00. The correct asymptotic form to use in this 
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region is 

I ig (X) = ^ZL (-e-^+ x + ie- x ), 

V^Fx V ) 

because the argument of X approaches — n from below 
(Abramowitz and Stegun 1972). In order to satisfy the 
boundary condition for large negative x, the coefficient of 
exp(— x) must be zero. But this condition does not lead to 
a constraint on x c . Instead, we obtain 

h/ ai = e" 2 ^, (24) 

which is incompatible with Eq. (|23[) . For large /3 the ratio 
(|24p is tiny. This very small ratio compensates for the se- 
vere damping that the left-travelling wave lip suffers in this 
region (see Eq. (pi"]) ). 

The requirement for two different linear combinations 
precludes the existence of a continuous analytic solution. 
There are no normal modes: when a 7^ the solution is not 
continuous at x = 0, when a = the solution is continuous, 
but possesses undefined derivatives at x = 0, a difficulty 
that can not be resolved whatever the sign of a. 

Let us attempt to understand this mathematical diffi- 
culty in more physical terms. The reason for the two incom- 
patible linear combinations is the heavy damping caused by 
the singularity near x = 0: this strong damping prevents 
the simultaneous satisfaction of the two boundary condi- 
tions. The correct combination of waves in the region x > 
violates the boundary condition in the x < region and vice 
versa. For example, suppose we take the correct combination 
of components in x > and analytically continue this so- 
lution into the x < region (in so doing passing above the 
singularity, as always). When our solution passes through 
the vicinity of corotation, the (right-going) I-ip component 
will abruptly pick up an 'amplification' factor (because 
we are proceeding 'backward in time' with respect to this 
wave), while the J„ component picks up a damping factor 
e~ w/3 (as we are proceeding forward in time). As a conse- 
quence, the amplitudes of the two waves will differ by e~ 2wl3 . 
But this large difference will upset the balance necessary to 
satisfy the decaying boundary condition at negative infin- 
ity, which requires the two components to be of comparable 
amplitude (at the very least). Thus as x becomes large and 
negative, u will blow up because it will be dominated by the 
coefficient of e~ x that issues from I v . If we take the correct 
combination of waves in x < and follow them into x > 
the same problem is encountered. 

For small (3 (less than 1/2), the order of the Bessel func- 
tions v is real and each solution component picks up a phase 
shift, rather than an amplitude change, when it is followed 
through corotation (see Eq. (|21[) ~). But this phase shift is 
also sufficient to prohibit the formation of normal modes as 
before: if we remove the divergent component of the solution 
when x — » 00, the phase jump across x = means we are 
unable to remove the divergent component of the solution 
at x — > —00. 

Though there exist no analytic eigenfunctions, when 
a — a piecewise solution can be devised which is con- 
tinuous at x = but whose derivative, as already noted, is 
undefined at this point: 

u B = Q(x)u+(x) + Q(-x)u-(x), (25) 

where Q(x) is the Heaviside step function and the right and 



left functions are 

u+ —alx 1 ^ 2 [h(x) — I- V (x)] , (26) 
u- =b{x 1/2 \l v (x) + e 2 ™I-v(x)\ , (27) 

where a\ and bl are two undetermined constants. In fact, 
because of translational symmetry, we can define such a 
neutral singular 'mode' at each x. As an ensemble this set 
of structures (understood as 'weak' or 'distributional' solu- 
tions) defines a continuous spectrum, because each radius 
is a singularity for at least one mode. Associated with the 
non-normal operators which appear regularly in shear flow 
problems, such 'solutions' aid in the evolution of initial data 
and thus only have physical meaning as kernels in the in- 
tegral expressions of u (see for example Schmid and Hen- 
ningson, 2001). In particular, the continuous spectrum is of- 
ten responsible for transient growth effects. The initial value 
problem we address in Appendix C, and in Section 5 where 
we show the relationship between the continuous spectrum 
and shearing waves. 

In summary, this simple demonstration shows the main 
problem trapped inertial waves suffer if they are to form 
around their corotation singularity. The shearing sheet 
model, though probably inadequate to fully describe iner- 
tial waves in real relativistic disks, here suffices to illustrate 
the principal physical and mathematical issues at play. In 
short, the absence of normal modes stems from the selective 
wave damping at the corotation radius. As a consequence, 
potential standing waves, composed of both right and left- 
travelling waves can never be calibrated so that they decay 
in the far field on either side of corotation. (In Appendix 
A these results are generalised to a slightly incompressible 
fluid, and in Appendix B the influence of viscosity is consid- 
ered.) 



4 ABSENCE OF TRAPPED INERTIAL WAVES 
IN A SEMI-GLOBAL MODEL 

In this section, we present a WKBJ calculation in a more 
realistic semi-global model of a disk. In the limit of large 
fe z , the main result of the previous section is recovered. The 
analysis is similar to, but somewhat simpler than, the de- 
velopment in K02. Nevertheless, it is still involved and the 
reader may skip directly to Section 5 without loss of conti- 
nuity. 

Our strategy is to first assume that trapped standing 
waves grow, a > 0, then to determine their approximate 
profiles near the two Lindblad resonances (subsections 4.2.1 
and 4.2.2 and 4.3) and near the corotation resonance (sub- 
section 4.2.3), and finally to match these three solutions in 
the WKBJ regions which border their domains of validity 
(subsection 4.4). The matching conditions supply an eigen- 
value equation (|58[1 for the growth rate a of the trapped 
modes. We approximately solve this equation, as in K02, 
and determine that in fact a < 0. We conclude, as a conse- 
quence, that no such trapped modes can exist. 

4.1 The wave-mode equation 

We begin with an incompressible disk in cylindrical geom- 
etry (r, cf>, z) orbiting with frequency profile f2(r). If the 
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isothermal approximation is used, the dynamics remain es- 
sentially unaltered (though more mathematically involved). 
The system admits an equilibrium state characterised by 
the (constant) density po and the flow u = rfle^. A small 
perturbation taking the form oc exp(im<f) + ik z z — iujt) is 
governed by the following linearised equations, 

iwu'r + 2S7ii0 = d r h' , (28) 

iutu'tj) — k 2 / (2S7) u' r = i(m/r)h' , (29) 

iuju'z = ik z h' , (30) 

= u r jr + d r u' r + i(m/r)u'^ + ik z u z , (31) 

where h = p' / po and u> — to — mil. The epicyclic frequency 
is defined through ft 2 (r) = 2Q,{rd r Q + 20) but, in order to 
approximately account for general relativistic effects near 
a black hole, «(r) may take other forms (see Okazaki et 
al. 1987, Kato 1990). 

This set can be manipulated into a simple second or- 
der ODE for u' r , which may be significantly simplified (as in 
K02) by assuming that the mode oscillates rapidly in com- 
parison with the radial variation of the background equilib- 
rium. If 

d r In k, d r In Q ~ K, 

where K is some characteristic wavenumber, then K <C k z . 
Also we let 

u r /r <C d r u r ~ k z u r , k z S> m/r. 

If we normalize the unit of length so that K — 1, the wave- 
shape equation is approximated by 



where 



d r u' r + e f(r, uj) u' r = 



(32) 



(33) 



and e = k z /K is a large dimensionless parameter. Equation 
((3"2"|) is analogous to Eq.@ in Section 2 and Eq. (20) in 
K02. For notational ease, hereafter we drop the prime and 
r subscript on u' r . 

Equation (J32J possesses two turning points _Ri and R2, 
the inner and outer Lindblad resonances, which are defined 
through 



respectively. These two radii circumscribe a resonant cav- 
ity in which low frequency standing waves may localise. In 
addition, there exists the corotation singularity R c defined 
through ui = 0. As noted by Kato (2001b, 2002), for nonzero 
m, this point usually falls between Ri and R2, though there 
exists a small interval of ui for which this is not the case (see 
Perez et al. 1997). This special circumstance is not investi- 
gated here and we assume Ri r < R cr < R^r (where the r 
subscript indicates real part). Because the resonant cavity 
contains the corotation singularity, neutral modes are pro- 
hibited. Consequently, if trapped normal modes are to exist, 
they must occur in growing/decaying pairs. 

We set uj = Lu r +ifj where uv and a are real and < a <C 
Q. The analysis is hence limited to slow growing modes. As 
a consequence, Ri, R2 and R c gain small imaginary parts. 



These can be approximated, to leading order in a, by 

- S5^« (34) 

R 2i = „° , , , (35) 



Rci — 



mfi' 



(36) 



where a prime indicates differentiation with respect to r, 
and a subscript of 1, 2, or c indicates evaluation at r = Ri r , 
7?2r, or R cr respectively. Even if k is assumed to possess 
a turning point near R cr (as it will in the inner part of a 
general relativistic disk) it is easy to see that Ru, Rn, and 
Rd are negative for all m and that 

Rli < Rci < R2i- 



4-1.1 'Rossby-wave term' 

Before continuing we draw attention to the 'Rossby-wave 
term' that will become important when some of the scaling 
assumptions are relaxed — namely, if we let k y ~ k z , and 
suppose there exist localised regions where the orbital fre- 
quency varies rapidly. If this were the case, the function / 
would pick up the following term, 



(m/r)d r (K 2 /n) 
2k\u 



(37) 



This simple pole may give rise to 'Rossby wave instability' in 
two-dimensional (k z — 0) dynamics (Lovelace et al. 1999, Li 
et al. 2000) , though for the large k z scalings we employ here 
Rossby modes do not appear. In particular, near corotation 
uj = the Rossby wave term is dominated by the double 
pole in (|33|l . Having said that, it is conceivable that when 
k z is small and the orbital frequency possesses a localised 
region of exceptionally strong shear, instability could arise. 
If so there will exist a critical k z above which the double 
pole dominates the dynamics and consequently eliminates 
potential Rossby wave instability. We do not attempt cal- 
culate this critical vertical wavenumber here, though we do 
investigate similar behaviour in Section 6 with respect to 
the Papaloizou-Pringle instability in a slender torus (see also 
Sternberg et al. 2008). In any case, realistic Q and k profiles 
are unlikely to exhibit such strong variation; if Rossby-wave 
instability is to occur it will probably have to rely on steep 
density gradients (which we have excluded with the assump- 
tion of incompressibility) . Q 



4.2 Approximate solutions near the Lindblad and 
corotation resonances 

4-2.1 Inner Lindblad radius 

Near r = R\ we can let / = ai(r — Ri) to leading or- 
der, where Qi = d r f evaluated at r = Ri. The wave- mode 
equation subsequently takes the form of a complex Airy's 
equation 



drU + e 2 «i (r — Ri ) u = 0, 



1 We thank the reviewer for bringing this issue to our attention. 
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with solution 

u = 01 Ai\J3i(Ri - r)] + h Bi[f3i{Ri - r)], (38) 

in which Ai[z] and Bi[z] are the two Airy functions, oi and 
6i are two complex constants, and 



a 2/3 / \-2/3 

Pi=e ' ai(-ai) 



(39) 



The constant ot\ will possess a small imaginary component 
whose sign plays an important role in negotiating the Stokes 
phenomenon we encounter later. For small a we can approx- 
imate the imaginary part of ai with the compact expression 



Im(ai) 



2a 



In 



■mfi' 



(40) 



As earlier, the subscript 1 indicates evaluation at r = Ri r . 
If the logarithm is expanded it is easy to see that the right 
hand side is positive. Then we can write 

a 2/3 1/3 27ri/3 

Because the imaginary part of a± is positive, the complex 
argument of /3i is a little greater than 2-7r/3 

4-2.2 Outer Lindblad radius 

Near r = R 2 , we can approximate the solution as in the pre- 
vious section. As before, we obtain a complex Airy's equa- 
tion 



d 2 u + 2 a 2 (r — R 2 ) u = 



(41) 



where a 2 = d r f, evaluated at r = R 2 . The solution takes 
the form; 



u = a 2 Ai[f3 2 (R 2 - r)] + b 2 Bi{[3 2 (R 2 - r)}, 

in which a 2 and b 2 are two complex constants, and 

a 2 /3 „ 1 „ n-2/3 

The imaginary part of a 2 to leading order is 



t r \ 2(7 ( 9 , 
Im(a 2 ) = — — m 

K2 V" r 



k' + mO,' 



(42) 
(43) 

(44) 



It can be shown that the right hand side is negative for all 
m for realistic k and f2 profiles. Now, 



02 



2/3, 



-012 



,1/3 



Because the real part of a 2 is negative, the complex argu- 
ment of (5 2 is a little less than — 7T. 



4-2.3 Corotation radius 

When r is close to R c we obtain the approximate equation 
d 2 u + k 2 a c (r — R c )~ 2 U = 0, 

with 



m 2 (f3') 2 

evaluated at r — R c . This admits a solution of the form 

u = a c (r - i? c ) 1/2+/5c + 6 C (r - Rcf' 2 '^ , (45) 

where a c and b c are complex constants and, to leading order 
in large e, 

/? c = e(-a c ) 1/2 . (46) 



We denote this solution hereafter by u c . 

The imaginary part of a c to leading order in small a is 
given by 



Im(a c 



Itdo 



d 



In 



(47) 



m 3 (ft' c ) 3 V* 

This can be demonstrated to be negative for all m > 0. Thus 

a -1/2 

4.3 Behaviour far from the resonant cavity 

We require the trapped mode to decay outside its resonant 
cavity, which is circumscribed by Ri and R 2 . By applying 
these two boundary conditions the two unknown constants 
b\ and b 2 can be determined. 

4-3.1 The region r < R\ r 

First we inspect the inner region of the disk, r < Ri r , far 
from the resonant cavity. In this region we have |e 2//3 (J?i — 
r)| ^> 1 for sufficiently large e, which supplies the following 
asymptotic forms for the Airy functions: 

AipiCRi - r)} = \Pi(R! ~ r)]" 1/4 (e- ( +ie^ 

Bi[/3i(7?! -r)] = [Pi(Ri -r)]- 1/4 (ie^ + e e ) , 



where 



e = ? [/3iCRi-' 



,3/2 



(48) 



Because of Stokes phenomenon the argument of /3i(Ri — r) 
is key and we must use (|40p to choose the correct asymptotic 
expressiorQ. This also tells us that Re(f) < 0. So in order to 
obtain a decaying it we must eliminate the exp(— £) terms 
and keep the exp(^) terms. This can be accomplished if 61 = 
i oi . The solution near the inner Lindblad resonance we now 
denote by iti and write as 



ui = oi { Ai[/?i(i?! -r)] + i Bi[/?i (i?i - r)] } . 



4.3.2 The region r > R 2r 



(49) 



We turn now to the outer region of the disk far from the 
resonant cavity and with r > R 2r . Here for sufficiently large 
e we have \e 2 ^ 3 (R 2 — r)\ 3> 1, and the Airy functions are 
approximated by 



Ai[(3 2 (R 2 - r )] = —= [/3 2 (R 2 - r)]~ 1/4 e" c 

2y/TT 

Bi[p 2 (R 2 -r)] = -j= [fc(R2 - r)]- 1/4 (l e~ c + e c 



with 



C=|[/32(i?2-r)f/ 2 . 



(50) 



The correct asymptotic form follows from Eq. (|44p . which 
also tells us that Re(£) > 0. To ensure decaying solutions 
we zero all the exp(£) terms, which corresponds simply to 



1 See for example http://functions.wolfram.com/Bessel-TypeFunctions/AiryAi 
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Radius 



Figure 1. A schematic diagram of the three solutions Mi, u c , and 
«2 in their regions of validity: near the inner Lindblad, corotation, 
and outer Lindblad radii respectively. These are represented by 
ill, Rc, and il2. Between the solutions are the regions W\ and 
W2 where the solution takes its WKBJ form. 

62 = 0. The solution near the outer Lindblad resonance is 
subsequently denoted by U2 and written as 

u 2 = a 2 Ai[/3 2 (i? 2 -r)]. (51) 



4.4 WKBJ matching 

In Fig. 1 the preceding three solutions are shown schemati- 
cally, each separated from its neighbour by either region Wi 
or W2- In these two regions we assume u takes its WKBJ 
form: 

u = aw Z" 1 ^ 4 exp e J / 1//2 dr^ 

+ bw f exp ^— ie J / 1//2 dr^ . 

Our task now is to extend solutions Ui, 112, and u c into re- 
gions Wi and W2 using the above ansatz and then to obtain 
matching conditions which set a. 

4-4-1 The iti solution in W\ 

Suppose that 7?i and R c are sufficiently spaced (or e suffi- 
ciently large) so that je 2/3 (i?i - r)| > 1 in region Wi- The 
correct asymptotic form for ui here is 

which can be manipulated into approximate WKBJ form. 
First we treat r as a complex variable, by simply adding to 



it a small imaginary component equal to iRn. Then we do 
the following: 

01 = r f 1/2 dr, 

= f Mr -El)] 172 dr + 0[(r - R 1 ) 5/2 ], 

J Rl 

^e" 1 + 0[(r - Ei) 572 ]. 

Here the path of integration is simply a straight horizontal 
line in the complex plane connecting Ri to r. We now have 

Ul « ai f~ 1/4 [e 1 ^ 1 + i e- 1 ^ 1 ) (52) 

in the region Wi, where ai is a new complex constant. 

4-4-2 The U2 solution in W2 

Suppose that |e 2//3 (i?2 — r)\ 2> 1 in region Wa- The correct 
asymptotic form for U2 here is 

^ = 1^2 MR* -r)]^ (e-<+^e<). 

As earlier, we treat r as a complex variable but this time 
add a small imaginary component equal to iJ? 2 i- Next we 
define 

02 = /' / 1/2 dr, 
JR 2 

where the integration contour is a straight horizontal line 
joining R2 to r. To leading order it can be demonstrated 
that e<f>2 — —i( and so 

U2 « 5 2 r lfi (e"" 02 + i e" 02 ) (53) 

in region VK 2 - 

4-4-3 The u c solution in Wi and W2 

So far our analysis has largely followed that of K02, but 
now we will deviate a little. The corotation radius introduces 
a problematic branch point in the asymptotic definition of 
u c and as a consequence it is obliged to take two different 
WKBJ expressions, one in region Wi and another in region 
W2- This subtlety appears to be incompletely treated in the 
two papers by Kato. Our approach here is relatively clunky 
but is unambiguous. 

The u c solution is rewritten as 

u c = Vr-^TTc (a c l ^~ R ^ + b c e~^ "*-*<>>) . 

Now to account for the solution in Wi we define 

03 = ^ f 1/2 dr, 

jR cr +iR u 

where the contour is the straight line connecting the point 
Rcr + i Ru to r in the region W\ . Here r has taken the small 
imaginary component i Ru . It can be shown to leading order 
in region W\ 

ie<f)3 = — p c hi(r — R c ) + eoi'J' 2 [iln(iAi) — n] , 

where Ai = Rd — Rn- Note the large damping/amplification 

1/2 

factor —ena c which appears to have been missed in K02. 
In addition, in Wi we have 

(r-i?,) 1 / 2 ^ 4 /" 174 . 
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Similarly in W2 we define 

^4= f f /2 dr. 

J R cr +iR 2i 

The contour of integration is the straight line connecting the 
point R C r + i R.2i to r in the region Here, however, r has 
taken a small imaginary component i i?2i ■ To leading order 

ie(f>4 — p c \n(r — R c ) — eaj i ln(— 1A2), 

with A 2 = -R C i — i?2i- For r in VK2, 



a/2 



«y 4 r l/4 . 



(54) 



(55) 



In summary, the WKBJ solution in region Wi is 
u c « i /~ 1//4 { 5 C exp (e-v/aT^ ln(iAi) — 7r] — *e0 3 ) 
+6 c exp (— e % /o7[iln(iAi) — 7r] + 16^3) j- 
The WKBJ solution in region VK2 is 

u c « ,f _1//4 { a c exp (iv / aleln(— iA2) + ec^) 
+6 c exp (— i v / aeln(— iA 2 ) — ie?^) j . 

4-4-4 Asymptotic matching 

First the two solutions (|52|l and (|54[) are compared in region 
Wi. We define 

nR cr +iR U 

#1 = / / 1/2 dr = 4> x - fa 

Jr.! 

and then equate coefficients of exp(iecii) and of exp(— iecii). 
This yields the two constants 2 C and 6 C in terms of 01 , which 
without loss of generality we let equal to 1. We have, as a 
consequence, 

a c — exp I en;!/ 2 [it — i ln(iAi)] — ie^\ j , (56) 
b c = -iexp |"eaj /2 [7r -iln(iAi)] + ie*i j . (57) 



Now we compare solutions (|53[) and (|55p in region Wi. 
We define 

*2 = ["* f 1/2 dr = 4> i -<t>2, 

J R cr +iR 2i 

and equate coefficients of exp(iec62) and exp(— it(f>2)- With 
Eqs (|56p and (|5T|) we obtains two equations, one which gives 
us 02 and finally the eigenvalue equation for uj. Once 5 2 has 
been eliminated and, after some algebraic manipulation, the 
eigenvalue equation takes the form 

'-Ai' 



1/2 



In 



A 2 



+ m 



(58) 



where n is some integer. In the next subsection we attempt 
to solve this equation approximately. 

4.5 Approximate solution to the eigenvalue 
equation 

For perhaps all realistic k and Q, the integral functions $1 
and 4 , 2 cannot be computed analytically. Instead, follow- 
ing K02 and K03a, we obtain approximations by using the 
asymptotic forms of / near Ri, R2, and R c . Doing so cer- 
tainly introduces some level of error, though the essential 
qualitative points we make should remain unaltered. 



The <&i integral is decomposed as follows 

r-Rcr+iR-li 

* 1 = / / 1/2 dr, 

Jr ± 

R 3 r R cr +iR u 

f 1/2 dr+ / f>*dr, 

Hi JR 3 

where R3 is an intermediate point on the line connecting Ri 
and R C i + iRu- In the first integral we approximate / by 
its asymptotic form near Ri, i.e. / = ai(r — Ri), while in 
the second we use its form near R c , f = a c (r — R c )~ 2 - The 
integrals are then straightforward to calculate and we find 

3 1 ' R 3r ~ Rlr) +Ci C In I J£ 



An analogous procedure furnishes us with 

2 
3 



* -^(-a2) 1/2 (^-i?4,.) 3/2 + ^ /2 ln(^^ 

-iA 2 



where we have introduced the intermediate point R4 which 
lies on the line connecting R2 and R cr + iR2i ■ 

With these expressions, the eigenvalue equation (|58[1 
becomes the more manageable 



■ 1/2 1/i, 
+ za„ 7r + a c In 



-a^Ai^ + ^(-a 2 ) 1/2 A.R 2 



i/2„ , „i/2,„ /ARsjMAi 
VAi? 4 -iA 2 



(59) 



where 



Ai?i = R 3r - Rlr, 
AR 3 = R cr - R 3r , 



AR 2 = 7?2r - i?4r, 
AR4 = i?4r — Rcr- 



Next we use the assumption that a is small and expand ai, 
a.2 , etc in powers of cr, 

ai = an + CK12 ia + 0(a 2 ), 

—a.2 = Q21 + «22 io + 0(a 2 ), 

Ctc = CKcl + «c2 icr + C((T 2 ), 

where the a coefficients on the right sides are all real and 
positive except for a C 2, (see Section 4.2). Note that K02 ne- 
glects these higher order corrections in the analogous terms 
of his analysis; however, these corrections will be impor- 
tant in the imaginary component of Eq. ((59} . We next take 
the imaginary part of Eq. (|59p but neglect the last term on 
the left side. This term contributes terms proportional to 
a]n(ARi/AR 2 ), a/ARi, and a/AR 2 all of which will be 
subdominant to the other terms if Ri, R2, and R c are suffi- 
ciently well separated so that the ARi are large and of the 
same order. Finally, we can solve for cr, 



o 1/2 



-1/2 A n 3/2 . -1/2 A „3/2' 

0120x1 Aity + Cl22«21 AR 2 ' 



(60) 



Thus a is negative. 

However, a negative a contradicts the initial assump- 
tion that a > (so crucial in negotiating the Stokes phe- 
nomenon). If, on the other, hand we begin the analysis with 
the assumption a < we finish with the opposite result: 
a must be positive. The only resolution to the impasse is 
to take (7 = 0. But this means that the modes will possess 
undefined derivatives. The conclusion is that there exist no 
trapped 3D normal modes for such WKBJ waves. 
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Before we continue, it should be noted that the negative 
results of Sections 3 and 4 and Appendix A pertain to spe- 
cific cases associated with certain approximations: e.g. the 
local shearing sheet model and incompressible and slightly 
incompressible fluids, and when in a more realistic geome- 
try, the WKB J limit of large k z . They are hence not proofs 
for the non-existence of all trapped 3D inertial waves near 
corotation. Together these results, however, offer a good ar- 
gument that this is indeed the case. 



5 SHEARING WAVES AND THE 
CONTINUOUS SPECTRUM 

In this section we briefly demonstrate the dynamical be- 
haviour of the continuous spectrum introduced in the shear- 
ing sheet analysis of Section 3. This is accomplished by in- 
verting the temporal Fourier transform, 



U(x,t) = J u(x,lo) 



' do;, 



(61) 



where u(x, u) is the solution to the modal problem and V is 
the appropriate integration contour in the complex uj plane. 
It is assumed that uj is real except near any singularity, at 
which point the integration path F will deviate above the 
singularity. 

For a typical initial value problem, in which the ini- 
tial condition is stipulated, we would construct u using a 
Greens function. Usually such problems throw up compli- 
cated integrals that are difficult to solve, and the analysis is 
usually limited to the large time asymptotic regime (Booker 
and Bretherton 1967, Watts et al. 2004, for example). In 
Appendix C we undertake such an analysis and show that 
localised initial conditions evolve towards decaying shearing 
inertial waves (for a description of shearing waves see John- 
son and Gammie 2005, Balbus and Hawley 2006). This long 
time behaviour issues from the collective influence of the 
corotation singularities. To emphasise this point, we present 
a slightly different calculation below in which each member 
of the continuous spectrum is summed equally. Doing so al- 
lows the inversion integral to be solved analytically and thus 
compared with the exact shearing wave solution of Johnson 
and Gammie (2005). 



5.1 Shearing waves 

Let us start with the integral (|61[) and immediately change 
the integration variable to 



fc_i_ x + 



k± uj 
qQ.k y 



(62) 



This is identical to the x* introduced earlier, but in the 
following it is denoted differently in order to avoid confusion 
with x. The new coordinate measures the distance between 
x and the corotation point specified by uj. In addition, a 
dimensionless time coordinate is introduced 



qQ.k y 



t 



We hence obtain 

/ qflky 



U = 



^ e %hAt)x f_u(x,0)e 



d9, 



(63) 



(64) 



where T is a contour in 6 space, u(x,uj) = u(x,6), and the 
time- dependent (shearing) radial wavenumber is 

k x {t) = qQ,k y t. 

A simple shift in time introduces a constant k' x into the 
above expression, and will render k x (t) in the familiar form 
used in Balbus (2003). The integral (|6ip now represents not 
a sum over all the uj but, for a given x location, a sum over 
all the corotation points that influence the dynamics. 

Before we specify the form of u the basic structure of 
the solution should be appreciated. First, we have a shearing 
wave contribution from the exp[ik x (t) x]. Second, this spatial 
dependence will be modulated by a possible ^-dependent 
factor issuing from the integral in (|64p. the details of which 
are linked to the initial condition. This integral will also set 
the time- dependence of the solutions' amplitude. For large 
times the amplitude decays like i -3 / 2± " (just as in Booker 
and Bretheton 1967). The dominant contributions to the 
integral in this case come from the combined criticial radii 
(the corotation points). 

To better clarify the action of the continuous spectrum, 
we treat an equal sum of its component singular modes; 
there is hence no modulation arising from localised initial 
conditions. Consequently, we set u(x, 0) = u s (9), where u s is 
defined in Eq. (|25|l . Now the spatial structure of the solution 
becomes a pure shearing wave with the integral in Eq. (|64|l 
depending on time only through r. We let this integral equal 
the amplitude function A(r), which we now calculate. The 
amplitude function involves a linear combination of the two 
integrals: 

f' oo 

At = / Vo [h(9) - /_„(<?)] e~ iT<> d0, 
Jo 

J — oo 



A 2 = 



dO. 



Both integrals are proportional to 



V6K V (6)< 



dO, 



which can be evaluated in closed form in terms of hyperge- 
ometric functions (Gradshteyn and Rhyzik 1963). Finally, 
after some tedious algebraic manipulation using a chain of 
functional identities, one can express the time-dependent 
amplitude of the shearing wave as 



A(T) = dF( 



' 2' 



2 ) 



3. 
' 2 ' " 



where ci and di are two constants and F(a, b, c; z) is a Gauss 
hypergeometric function, often represented as 2F1. This am- 
plitude is exactly that of a pure shearing wave computed di- 
rectly from the governing equation §8§ rendered in shearing 
coordinates: 



d 2 T A + ^{I + t 2 )- 1 A = 0, 



(65) 



(for the derivation see Johnson and Gammie 2005 or Balbus 
and Hawley 2006). 

The main point to take away is that inertial shearing 
waves are equivalent to the continuous spectrum issuing 
from the corotation singularities. This simple relationship 
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has not been generally appreciated in the astrophysics com- 
munity, but is quite natural (Craik and Criminale, 1986). A 
shearing wave understood this way is but a sum of structures 
Us each tightly localised to its critical radius and convected 
perfectly upon the background shear, like tracers. 



6 INSTABILITY IN A SLENDER TORUS 

The previous sections have investigated the role of the coro- 
tation point in an extended disk. The analysis was local 
and the disk-boundaries were sent to infinity. Consequently, 
trapped modes could not arise. In this section a different 
scenario is attacked: an inviscid slender torus, the radial 
boundaries of which play the essential part. Specifically, we 
examine the effect the corotation point has upon the in- 
compressible three-dimensional Papaloizou-Pringle instabil- 
ity (Papaloizou and Pringle 1984, 1985, Drury 1985, Goldre- 
ich et al. 1986). 

The classical model of this global hydrodynamical in- 
stability is two-dimensional, in which case the corotation 
point is not a singularity of the governing equations. The 
mechanism of instability proceeds in this case from the effi- 
cient communication of angular momentum from the inner 
boundary to the outer via wave motions across the corota- 
tion point. If, however, the mode in question posesses verti- 
cal structure, the corotation point, being singular, will im- 
pede this transfer and either inhibit growth or prohibit the 
formation of normal modes entirely (as in the previous sec- 
tions) . In the following we characterise this phenomenon for 
general q. Note that our results complement those obtained 
by Sternberg et al. (2008) who examine a similar problem 
with fixed boundary walls. 



6.1 Mathematical set-up 

We consider the shearing sheet approximation to a narrow 
incompressible fluid annulus with no equilibrium vertical 
structure. The centre of the sheet is anchored at the cen- 
tre of the annulus ro so that the fluid occupies the radial 
region circumscribed by \x\ ^ s, where 2s is the total radial 
width of the annulus. The Eulerian perturbation equations 
are Eqs (O-©, and we draw our general solution for u' x for 
a given (k y ,k z ,ui) mode from Eq. 

The appropriate boundary conditions for a narrow an- 
nulus is that the Lagrangian enthalpy perturbation must 
vanish at the free surfaces x — ±s. This condition may be 
represented by 



where 



,/ . a dho 

h + Ax —— = 0, 

ax 



Ax ■■ 



(66) 



qQkyX + ia 



is the radial displacement at the boundary. The Eulerian 
enthalpy perturbation is supplied by 

k u 



no 

and the radial gradient of the equilibrium enthalpy ho is 



This enthalpy gradient aids gravity negate the centrifugal 
force when the rotation law is non-Keplerian. 

In this section we consider modes that are corotating 
with the centre of the annulus, so that uo r = 0, but which 
possess nonzero a. We also scale space by fcj as in Section 
2.3, while introducing the scaled growth rate 

-- J^L ( — ) 

qk v \QJ ' 

(which in Section 2.3 is equal to — x c ). The boundary con- 
dition can now be expressed neatly as 

g 2 (a, + isf + {g(2 -q)(x + ia) 
ax 

-(2q-3)(l + ki/kl)x}u' x =0, (67) 
which must hold at x = ±s. The expression for u' x is 



u' x — A\f x + ia l v (x + ia) + By 'x + ia I- V (x + ia), 

where A and B axe constants to be determined, and u can 
be computed from (|16|) . Once A and B are eliminated from 
the two equations we arrive upon the dispersion relation for 
<j, which is a function of the three dimensionless parameters, 
q, k z /ky, and s. 



6.2 Growth rates 

Unfortunately, the dispersion relation is transcendental in 
a. In the interest of readability, the full expression is ban- 
ished to Appendix D. In general, the growth rates must be 
computed numerically. Nevertheless, some analytic progress 
can be made if it is assumed that the width of the annulus is 
small relative to the typical lengthscale of a (k y ,k z ) mode, 
i.e. |s| <C 1. If the annulus was much wider, especially if the 
mode's Lindblad resonances fitted into the domain, com- 
munication between the two boundaries via incompressible 
motions would fail and the instability would evaporate. In 
addition, we let the growth rates scale as a ~ s (Balbus, 
2003). This permits the scaling I v (ia±s) ~ s v , which helps 
us fish out the most significant terms in the dispersion rela- 
tion. To leading order it can then be boiled down to 



2q 2 (2~q)f(a 1 )aj + 2q 2 i(2q- 



3) vg(ai) a i 

3(3-4 2 ) + S(3-2<z) 
S 



0, (68) 



where we have set a = sa\ + 0(s 2 ) and where 

/fa) = (tffl - l) 2 ' y - (ttTi + if. 

g((Ti) = (ia! - l) 2 " + (iai + l) 2 " . 

Though the reduced dispersion relation cannot be solved 
analytically for a\ , it yields two significant results. 

First, we can establish the two-dimensional stability cri- 
terion easily. By setting k z — 0, and consequently v = 1/2, 
we get the twinned growing and decaying modes 



^ = ~(2q~3)tfx. 
ax 



qai = ± N /3(g 2 - 3) 

in agreement with Papaloizou and Pringle (1985) (see also 
Balbus 2003). The instability criterion is thence q > 

Second, it yields the curve of marginal stability for gen- 
eral three-dimensional modes once we set ai = 0. The curve 
is described by 

3 



Qcrit 



A{k z /k y ) 



(2(k Z /ky) 2 - y/3-{k Z /ky 



(69) 
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Figure 2. The marginal stability curve in the parameter space 
(k z /k y ,q) given by H69I I. The region above it is unstable. The 
region bounded by it and the vertical k z /k y = axis circumscribe 
parameters for which no modes exist. 

The curve of critical q as a function of k z /k y is plotted in 
Fig. 2. There we can see clearly that the greater the vertical 
pitch of the wavevector, the greater the critical q necessary 
for instability. 

As expected, the most unstable modes are two- 
dimensional. Three-dimensional modes, on the other hand, 
must contend with the singularity at corotation x = 0, 
which interferes with the mechanism of instability. As dis- 
cussed in Section 2.4, waves with small wavevector pitch, 
k z /ky ^ gfi/(2fc) (i.e. /3 ^ 1/2), pick up a phase shift upon 
crossing this point. But the singularity is far more destruc- 
tive for modes with greater pitch, k z /k y > gfi/(2/t) (i.e. 
(3 > 1/2) which are severely damped. No normal modes 
(growing or otherwise) of the latter class are possible in the 
slender torus. This is also true for some low k z /k y modes 
which only suffer the phase shift. Generally, however, such 
P ^ 1/2 modes do exist and occur as growing/decaying 
pairs. They are situated in the region above the marginal 
curve in Fig. 2. Their growth rates are lower than the equiv- 
alent 2D mode at the same q because the phase shift renders 
the transfer of angular momentum across corotation less in- 
efficient. This is illustrated particularly well in the eigen- 
function profiles in the next subsection. 

6.3 Eigenfunctions 

To examine the effect of the corotation singularity on the 
structure of the unstable modes we fix q = 1.8 and s = 0.1 
and subsequently vary the ratio k z /k y from to the value 
at which modes cease to exist (near 1.5). The eigenprob- 
lem Eq. (|D1|) is solved numerically by a Newtown-Raphson 
method. 



Figure 3. The modulus and argument of the u' x eigenfunction's 
radial structure for q = 1.8, s = 0.1, and various k z /k y . Panels (a) 
and (b) show the (two-dimensional) cigenfunction when k z /k y = 
0. The growth rate here is a = 7.60 X 10~ 2 Q. Panels (c) and 
(d) show the case when k z /k y = 0.5. Note that the presence 
of the corotation singularity breaks the radial symmetry. As a 
consequence, the unsealed growth rate is less: we compute a = 
6.23 X 10 — 2 f2. Finally, panels (e) and (f) show an eigenfunction 
with k z /k y = 1.4, a value which positions the mode very near 
the stability boundary sketched in Fig. 3. The unsealed growth 
rate is a = 9.91 X 10 — 4 Q. Note the near discontinuity in the 
argument at corotation. Motions on either side of this point are 
out of phase by a quantity near tc/2, and thus the transfer of 
angular momentum and energy from the inner boundary to the 
outer is severely impeded. The eigenfunctions are normalised so 
that the maximum of \u' x \ is 1. 



In Fig. 3, the top two panels (a) and (b) exhibit the 
Eulerian radial velocity component (u' x ) of the classical two- 
dimensional mode (k z — 0) as a function of dimensionless x. 
Both the modulus and argument are presented. The most 
important features here are the relatively small variation in 
radius (x), which is of order s 2 , and particularly the small 
change in phase from the inner edge to the outer edge. 

In the next two panels, Figs 3c and 3d, a three- 
dimensional mode is plotted with k z /k y — 0.5. The mode 
grows slower than the previous case, and the symmetry of its 
structure is broken. In Figs 3e and 3f the k z /k y = 1.4 case is 
presented, which is near criticality. Any value of k z /k y much 
larger does not return a modal solution. The salient feature 
here is the near discontinuity in Arg(u^) at the corotation 
radius. Motions on either side of this point are out of phase 
by a quantity that approaches 7r/2. This represents the limit 
of severe damping of energy and angular momentum transfer 
and thus zero growth. 
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7 CONCLUSION 

Axisymmetric inertial waves can become trapped in the in- 
ner regions of disks orbitting black holes, forming standing 
modes (Okazaki et al. 1987). In certain cases a small set 
of non-axisymmetric inertial waves can do the same (Perez 
ct al. 1997). However, most 3D non-axisymmetric inertial 
waves possess a corotation singularity within their resonant 
cavity and, as a consequence, these waves will have difficulty 
forming trapped standing motions. 

This idea was demonstrated with a simple local model, 
the shearing sheet, which crystallises the important physi- 
cal and mathematical points thus allowing a straightfoward 
interpretation of the problem. Its conclusions are bolstered 
by a WKBJ analysis in a more realistic semi-global model in 
cylindrical geometry in Section 4, and a compressible analy- 
sis in Appendix A. The central point is that a trapped stand- 
ing wave must carefully balance its two component travelling 
waves in order to satisfy the two decaying boundary condi- 
tions far away from the corotation region; at the very least, 
these two travelling waves must possess comparable ampli- 
tudes. But the corotation singularity acts as a powerful wave 
absorber, heavily damping one of the waves relative to the 
other and hence destroying the necessary balance. If one 
boundary condition is satisfied the other must be violated. 

In summary, non-axisymmetric 3D inertial waves should 
play little direct role in QPO models based on diskoseismol- 
ogy, though they may be implicated in parametric insta- 
bilities, where they can act as 'intermediaries' transferring 
energy between two trapped axisymmetric waves (for exam- 
ple, Ferreira and Ogilvie 2008). More generally, the viability 
of diskoseismology rests on the interactions between trapped 
inertial waves and their turbulent environment. Recent lo- 
cal and global simulations, in particular, have shown that 
trapped modes struggle to emerge from MRI-induced tur- 
bulence (Arras et al. 2006, Reynolds and Miller 2009). 

The corotation singularity, while forbidding discrete 
normal modes, on the other hand generates a continuous 
set of neutral singular 'modes', a continuous spectrum. We 
show that an ensemble of these singular structures corre- 
sponds exactly to the inertial shearing waves computed by 
Johnson and Gammie (2005) and Balbus and Hawley (2006), 
and in Appendix C we use them to solve the initial value 
problem for localised initial conditions in the asymptotic 
limit of large time. This provides an alternative, and an- 
alytically tractable, interpretation which helps deepen our 
understanding of shearing waves, which appear naturally in 
shearing box simulations and may be important in planet- 
disk interactions (Balbus 2003, Balbus and Hawley 2006, 
Shen et al. 2006). 

Lastly, unstable 3D modes in a simple model of a slen- 
der torus were studied in the context of corotational damp- 
ing. As expected, 3D modes of small pitch k z /k y grow at 
slower rate than their 2D counterparts on account of the 
impeded angular momentum transfer at the corotation ra- 
dius, which is a singularity when k z 7^ 0. Modes of greater 
pitch do not exist at all. This simple example analysis fleshes 
out rather nicely the ideas put forward in Li et al. (2003). 
Moreover, it shows that astrophysical tori which are unsta- 
ble to the incompressible Papaloizou and Pringle instability 
will be dominated by the two-dimensional unstable modes 
above all others. 
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APPENDIX A: DEPARTURES FROM 
INCOMPRESSIBLITY 

This appendix reinforces the result of Section 3 by inves- 
tigating the role of compressibility. The most salient effect 
of compressibility in our problem is the 'leakage' of iner- 
tial waves through the walls of their confining potential 
wells. We establish whether this tunneling can circumvent 
the strong constraints imposed by the boundary conditions 
and consequently permit trapped normal mode solutions. 
Recently Ferreira and Ogilvie (2008) numerically demon- 
strated that a trapped axisymmetric inertial mode can tun- 
nel through its confining barrier and emerge on its other 
side as a small-amplitude p-mode. The effect seems small 
but may be important in the normal mode formation (or, 
rather, non- formation). 

In light of this we undertake a generalisation of the anal- 
ysis of Section 3 to a 'slightly incompressible' fluid; that 
is, we define a small parameter associated with the sound 
speed and expand to its first order. The incompressible anal- 
ysis of the main paper is subsequently interpreted as the 
zeroth order approximation. As before, we attempt to de- 
scribe trapped normal r-modes but again we find that no 
such solutions exist. 

Al Governing equations 

We return to the set of linearised equations for a small dis- 
turbance but now use the full continuity equation 

dtp + u • Vp = — pV • u, 

instead of the incompressibility restriction V ■ u = and 
we assume an ideal fluid, so that p = c 2 s p, where c s is the 
(constant) sound speed. In addition, the tidal potential $ 
picks up a component representing vertical gravity. 

The vertical equilibrium of the basic state exhibits 
a Gaussian vertical structure. Consequently, we assume a 
mode structure proportional to 

exp(ifc y y — iuit) Re„[z/(V2H)], 

where He„ is a Hermite polynomial of order n and H — c s /Q, 
is the disk scale height. In the shearing sheet model this 
ansatz yields the linearised equations, 

itch! = c? a [d x u x + ikyUy + i(k 2 /u>) h'], (Al) 
iuiu' x = ~2Q,u y + d x h', (A2) 
iuiu'y = ^(k 2 /Sl)u' x + ik y h' , (A3) 

in which we have defined k z = i/n/H, and eliminated u' z . 
This set of equations can be combined to form an equation 
for u' x , 

2 / 2q£lujk y i 

(TKj_ — CO 



This equation exhibits both the regular singular point at 
corotation, lo — 0, and the apparent singularities at the ver- 
tical resonances, uj 2 — (? a k\. The latter points, with the 
Lindbland resonances uj 2 = k 2 , circumscribe the two poten- 
tial barriers. We set u) = uj t + i<? and change coordinates 
as in Section 2. Doing so introduces the new dimensionless 
parameter 

8 = (Hk±y\ 

which measures the importance of compressible effects. To 
make further progress the key step is to assume this quantity 
is small and to expand the mode equation in its powers. The 
result up to order S 2 is 

d 2 u + 28 2 a 2 XdxU 

+ {|i-(l-7)+a 2 A 2 j M = 0, (A4) 

where X = x — ix c , with f3 and x c given in Section 2, and 

a = q(k y /k±)S, 7 = 8 2 K 2 [q(k y /k ± ) 2 - 1]. 

Note that we have dropped the prime and subscript x on 
u'x- In addition, we could expand u in small 8 as well, but 
we retain all its subdominant components. 

The equation ()A4[) is governed by the dimensionless 
parameters q, k y /k z , and 8. By setting 8 — we recover 
Bessel's equation (|12p . The additional terms alter the nature 
of the potential, which is described by the negative of the 
term in curly brackets. To leading order, the two forbidden 
zones are contained, on the one hand, between x = f3 and 
x — oT 1 , and, on the other, between x = — /3 and x = —a" 1 . 
In the incompressible limit we have 5^0 and find the bar- 
riers extend to infinity, as in Sections 2 and 3. There are thus 
three zones in which wave-motion can take place. Modes lo- 
calised primarily to the two regions \x\ > a -1 , far away 
from corotation, we identify as p-modes (density waves), 
and modes localised near corotation \x\ < [3, we identify 
as r-modes (inertial waves). 

A2 Mathematical analysis 

Equation ()A4|) can be solved analytically by assuming the 
functional form 

u = X a e bx2 f(X), 

where a and b are constants we are free to choose. Af- 
ter substitution of this ansatz into ()A4|) a and b are set 
to values which zero the x~ 2 and x 2 terms in the coeffi- 
cent of /. The result is a version of Rummer's equation 
which can be made more explicit by the coordinate trans- 
form Y = iay/1 — a 2 x 2 . Then we have 

Y&f + (l±v-Y)dYf- V f = 

where v was introduced earlier in (|16[1 and 77 is given through 

_ l±v i(l + a 2 - 7) 

^ 2 4(3^1 — « 2 

Note the ambiguity in the sign of v which arises from a free- 
dom in the choice of a and which corresponds physically to 
the direction of the p-mode solution. Without loss of gener- 
ality we take the positive sign (the left-going p-mode), but 
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Figure Al. Here we present the real part of u± for x > with 
parameters q = 3/2, k z /k y = 7, and 8 = 0.425. In addition, we 
indicate the location of the resonances: the corotation point is at 
x = 0; the solid circle refers to the Lindblad resonance, x = 4.67, 
and the asterisk refers to the vertical resonance, x = 11.1. The 
forbidden zone is hence delimited by an asterisk and a circle. We 
only plot Mi to identify it as possessing the features of an inertial 
wave. As is made clear in the text, inertial waves do not manifest 
as normal modes in this problem. 

remain conscious that results pertaining to it also hold for 
the negative case (right-going mode). 

The two independent solutions to Kummer's equations 
are M(rj,v + 1,Y) and U{rj,v + 1,Y), the confluent hy- 
pergeometric functions of first and second kind respectively 
(Abramowitz and Stegun, 1972). The two linearly indepen- 
dent solutions to the original wave-shape equation (IA4|) . 
with the correct values of a, b, and Y substituted, are hence 

Ul = X 1/2+ " exp (-Iq 2 X 2 - | i a^T^X 2 ) 

x U(r), l + v; iaVl- a 1 X 2 ) (A5) 

u 2 = X 1/2+ " exp (-±a 2 A 2 - § i a^Jl-cPX 2 ^ 

x M(rj, l + is; ia^/l- a 2 X 2 ). (A6) 

We identify solution Ui with the r-mode and ui with the 
(left-going) p-mode by inspection of their profiles (see Figs 
Al and A2): mi is confined almost fully to the cavity near 
corotation with some leakage through the potential barriers; 
U2 is confined to the areas outside the potential barriers with 
some leakage into the corotation region (which is too small 
to see). Both solutions decay speedily for large \x\. 

A3 Absence of trapped inertial waves 

First, it is assumed that a > 0. Next, consider the argument 
of the hypergeometric function in ui , 

Y = iay/1 — aX 2 = u\J 1 — a 2 \~2ax + i(x 2 — a 2 )] . 

Suppose that we following the solution along x towards the 
corotation point from the right (in region x > 0). When 
we reach x = a the imaginary part of Y changes sign but 
the real part of Y remains negative, which means that in 



Figure A2. In this figure we plot the real part of the left- 
travelling p-mode as described by U2 as a function of x. The 
parameters are q = 3/2, k z /k y = 1.75, and 8 = 0.1. The Lind- 
blad resonances are at x = ±1.17 and are denoted by solid circles, 
the vertical resonances at x = ±13.4 and denoted by asterisks, 
and the corotation point is at x = 0. 



the complex K-plane we cross the negative real axis. But 
the negative real axis happens to be a branch cut for the 
function U (??, v+l,Y) (Abramowitz and Stegun 1972). This 
means that 

lim U(ri, v + 1 ,—Z + iS) ± lim U(rj, v + 1, -Z - iS) 

6^0+ 6^0+ 

for any Z positive and real. So, in effect, when a 7^ the 
solution ui exhibits a discontinuity when x = a and is not 
analytic as a result. (In fact, wi proceeds to a different Rie- 
mann surface.) If we set a = then the discontinuity is re- 
moved but then the solution possesses undefined derivatives 
at x = and is thus not analytic either. Physically, the Ui is 
a combination of left and right-travelling p-modes, organised 
so that |«i| — > as jx| — ► 00. But just as in Section 3, the 
satisfaction of these two boundary conditions comes at the 
price of continuity — of the function itself or its derivatives. 
In summary, there can be no trapped inertial modes even 
when compressibility is taken into account. 



A4 Density waves 

On the other hand, the confluent hypergeometric M of the 
p-mode U2 is an entire function and thus suffers no branch 
cuts. If there is a small nonzero a the p-mode is analytic 
throughout the domain: from Fig. A2 it impacts on the po- 
tential barrier from the right, decays exponentially in the 
forbidden region, and trickles into the narrow corotation re- 
gion before emerging as a smaller amplitude p-mode oscilla- 
tion on the far side of the left-most barrier. We suspect the 
numerical integration undertaken in the Appendix A of Li 
et al. (2003) produces a profile not dissimilar to uz, though 
it may not capture the function's far-field decay. 
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APPENDIX B: INFLUENCE OF VISCOSITY 

Typically the continuous spectrum gives rise to algebraic 
growth or decay. However, even a small amount of dis- 
sipation will remove the singularity at corotation, and in 
bounded flows the inviscid continuous spectrum will congeal 
into a discrete but dense set of eigenvalues each associated 
with a regular eigenfunction exhibiting an exponential time 
dependence. In the unbounded shearing sheet, however, the 
analogous spectrum will remain continuous and its eigen- 
functions singular. 

Though the singularity at corotation is removed by vis- 
cosity, an incident travelling wave should still suffer strong 
absorption when it strikes the corotation region. This means 
the arguments of Section 3, which we used to dismiss the 
possibility that left and right-travelling waves can be com- 
bined to form trapped standing waves, still apply qualita- 
tively. On the other hand, one can demonstrate the im- 
possibility of constructing a localised time-exponential so- 
lution from a sum of viscous shearing waves, and thus of 
viscous standing eigenmodes generally (Gordon Ogilvie, pri- 
vate communication). It follows that the unbounded shear- 
ing sheet must possess a viscous continuous spectrum which 
describes the dynamics of viscous shearing waves. 

Just as in the Orr-Somerfeld problem, the effect of a 
small viscosity in a disk is most pronounced at the critical 
radius - at corotation. A viscous boundary layer develops 
around this point to mediate the absorption of wave energy. 
From the viscous equations of motion we can derive the fol- 
lowing sixth-order equation in the dimensionless variables of 
Section 2.3, 

E 2 d 2 u- E 2 u + 2 u = (Bl) 
where the operator E is defined by 

S= (x-ix c )-iTZ' 1 (d 2 -l), (B2) 
which introduces the (modified) Reynolds number 

*-(£)&)■ 

with v the kinematic viscosity. From this equation it can 
be shown that the internal viscous layer at corotation pos- 
sesses a width of order TZ -1 ^ 3 when the Reynolds number is 
large. Approximate solutions very close to corotation (inside 
the viscous layer) can be easily derived for WKBJ waves 
when f3 ~ TZ 3> 1. All six solution decay (or grow) expo- 
nentially through the layer at a (spatial) rate (1Z /3) 1//3 . This 
means the total reduction in amplitude of an incident WKBJ 
wave through the entire layer will be of order /3 1//3 , which 
should be contrasted with the rate /3 that emerges from the 
inviscid analysis of Section 2.4. Thus, viscosity allows in- 
cident waves to penetrate corotation with greater success, 
though the damping will still be sufficient to preclude stand- 
ing eigenmodes. 



APPENDIX C: THE INITIAL VALUE 
PROBLEM 

This appendix offers a more formal treatment of the initial 
value problem than in Section 5. We return to the original 
linearised equations ©-([Til and take Fourier transforms in 



space and time. The temporal Fourier transform is defined 
by 

roc 

f= / fe^dt. 
Jo 

It follows that 

— iuju x — 2Q,u y = —d x h + u x , 

- K 2 - ., f 

— lUJUy + — Ux = —IKyh + U y , 

— iuJUz = —ik z h + u a z 

with the incompressibilty condition. It is assumed that the 
initial conditions, u° and h°, are analytic and obey the far 
field decaying boundary conditions. Like in Section 3 this 
set is reduced to a single equation for u x . With the space 
dimension scaled by fe^ 1 and with the transformed velocity 
field scaled by fi, we obtain 

^ + {j£^- 1 )^ = jzhj A{x > x ° ) > (C1) 

where x c = — (ujk±_) / (qQk y ) . Note that this definition is dif- 
ferent to that appearing in Sections 2 and 3. Equation (|C1[I 
is the inhomogeneous Bessel's equation with 

A(x, x c ) = i k ^ (d 2 u° x - u°) + . ru°. 

QKy Q [X x c ) 

To simplify the following mathematics, without altering our 
conclusions, u° z is set to 0. Thus A = A(x). 

CI Greens function solution 

Equation (|C1|I can be solved using a Greens function, 
G(x, y), so that 

J -co y-xc 

in which the integration contour deviates below the singular- 
ity at y = x c - The Greens function must satisfy the equation 

where 5{x) is the Dirac delta function, while at the same 
time satisfying the far field decaying boundary conditions. 
A suitable choice is 

G(x,y) = —{Q(y-x)-u+(y)u-(x) 

+Q(x-y) -u+(x)u-(y)}, 
where Q(x) is the Heaviside step function, u+ and it- are 

u + — (x — x c ) 1/2 K v {x — x c ), 
u- = (x— x c ) 1/2 (I v (x — x c ) — e 2 """I-„(x — x c )) 
and W is their Wronskian: 

We can now write down the solution to Eq. (|CT|| , 

Ux=Au- + Bu+, (C2) 
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where 



1 f c 

B(x,x c ) = ± J 



h{y)u + (y,x c 
y-Xc 



dy, (C3) 



1 f x k{y)u-{y,x c ) 
y - %c 



dy. (C4) 



C2 Evolution at large times 



Expression (IC2[) is put in the inversion integral to achieve 
the full time-dependent solution of u' x . From (|64|l in Section 
5 this can be written as 



U = 



q£lk y 
~k~ 



ik x (t)z 



u(x, 9) ( 



d6, 



(C5) 



where 9 = x — x c and F is a suitable contour in 9 space. 
Needless to say, general initial conditions do not yield closed 
forms for A and B, let alone U. However, some progress can 
be made if we take the asymptotic limit of large time t> 1. 
According to the Riemann-Lebesgue lemma the dominant 
contribution to the 9 integral occurs when 9 is small because 
of the exp(— ir9) factor. Consequently, we expand A, B, u+, 
and it- in 9. For A and B this means x c ~ x, and so A = 
A(x) and B = B(x). If £ = y — x is a dummy variable, we 
get 

A{x)* w J A( I + Qr 1/2 [i40-^KM. (C6) 

and a similar expression for B. We assume for the moment 
that these leading order expressions for A and B are nonzero 
and well-defined. 

With A and B no longer depending on x c (and hence 9), 
we can tackle the 9 integral just as in Booker and Bretherton 
(1967). For small 9 this integral is proportional to 

I (c{x)9 1/2+l/ + D(x)9 1/2 -^ e' lT °d9, 

where F e is that portion of the integration contour near 9 — 
and C(x) and D(x) are functions involving combinations 
of A and B. An order 1 variable is introduced, £ = ir9 which 
establishes the basic time dependence of the integral. It can 
as a consequence be re-expressed as 



C(x)r 



-3/2+1/ 



+ D(x)t 



-3/2-1/ 



where the new functions C and D involve integrals of the 
form 



c 



l/2±u -C 



e~ k dC, 



which may be integrated numerically. 

Finally, we write down the long time fate of a localised 
initial condition. Given k y and k z we have 

u' x oc e ik ' ll),rt " ,+ ' 1 '' t~ 3/2 (C(x) t u + D{x) r") . (C7) 

The mode decays algebraically, with oscillatory behaviour 
when v is imaginary (i.e. (3 > 1/2). The spatial structure of 
the solution is a shearing wave localised within an 'envelope' 
defined by the functions C and D, which in turn depend on 
the initial condition selected. 



C3 The behaviour of the A and B functions 

To derive the main result we assumed that the leading order 
terms of A and B are nonzero and well behaved when 9 is 
small. In this subsection we attempt to justify these assump- 
tions. In (|C6fl we expand A(x + £) in a Taylor series around 
x and truncate at some order N. This supplies a reasonable 
approximation as the dominant contribution to the integral 
comes from near £ = 0; moreover u+ decays rapidly with £ 
(like e - ^). The first term in the expansion is computed, and 
we find 

Ajx)^^^^) 

A similar expression exists for B. Thus the leading order 
terms exist and are regular. We can thus be assured that 
generally A and B axe both nonzero and well behaved. 



APPENDIX D: DISPERSION RELATION FOR 
INSTABILITY IN A THREE-DIMENSIONAL 
SLENDER TORUS 

The full dimensionless dispersion relation which issues from 
the solvability of (|67[) at x = ±s can be expressed as 

X4(ff) • + X3(ct) • ct 3 + X2(ct) • a 2 

+ Xi(v)-v + Xo(°)=0, (Dl) 

where 

X4 = 4<j 4 {S\+l - Situ) , 

X3 = 2iq 3 { (q-4- 2»q){Sl-» - S^ v ) 

+ (q-4 + 2v q )(S^-Sltn}, 
X 2=8q 2 (2-q)(l + k)(Sl u -S- v ) 

+ 2sq 2 (Q -8q + q 2 + fe(6 - 4q) - 2uq 2 )(Sl-" + S?_„) 

- 2sq 2 (Q -Sq + q 2 + fc(6 - 4q) + 2vq 2 ){S 1 +J J + S^+ u ) 

+ 8q 4 s 2 (Sl+" - SltZ), 

X i = -2iq 2 s { 4(2q - 3)(1 + k)v(Sl v + S'") 

+ s(12 - 4q - q 2 - 4k(2q - 3) + 2vq 2 ){S 1 - v - S v x _ v ) 
+ s(12 - 4q - q 2 - 4k(2q - 3) - 2vq 2 )(S 1 + l/ " - S£ v ) }, 

Xo = s 2 {4[9 - 3q 2 +k(3- 2q) 2 ](l + %)(Sl„ - SZ") 

+ 2q 2 s{-6 + q 2 + k(4q - 6) - 2uq 2 ) + S 1 '") 

- 2q 2 s(-6 + q 2 + k(4q - 6) + 2uq 2 )(S^„ + S l _+ V ) 

-ti-i 
'l+i- 



in which k — (k z /k y ) 2 and where we have used the short- 
hand 

Sv = I + s) I v (ia — s). 



When a ~ s <C 1, Equation ()D1|I reduces to (|68p after using 
the scaling 



Consequently, we find only the last three terms in (IDlj) con- 
tribute to the leading order balance. 



